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A GENERAL ALGORITHM USING FINITE ELEMENT METHOD FOR 
AERODYNAMIC CONFIGURATIONS AT LOW SPEEDS 

By 

R. Balasubramanian^ 

SUMMARY 

A finite element algor ithr for numerical simulation of two- 
dimensional, incompressible, viscous nows has been developed. The 
Navier-Stokes equations are suitably modelled to facilitate direct 
solution for the essential flow parameters. A leap-frog time 
differencing and Galerkin minimization of these model equations 
yields the finite element algorithm. The finite elements are 
triangular with bicubic shape functions approximating the solution 
space. The finite element matrices are unsymmetrically banded to 
facilitate savings in storage. An unsymmetric L-U decomposition is 
performed on the finite element matrices to obtain the solution for 
the boundary value problem. 


INTRODUCTION 

Under grant NSG-1094, a computer program for mamerical solution 
of two-dimensional, incompressible, viscous flows has been developed. 
The program is operational and has been used to simulate the "driven 
cavity" problem. The numerical results for the "driven cavity" are 
compared with known finite difference solutions. The purpose of 
this report is to provide documentation of the preliminary numerical 
results as well as to present the software programs generated under 
grant NSG-1094. 

The Analytical Fluid Mechanics Section, FMB, HSAD, LaRC has been 
in the process of evaluating the efficiency of finite element 
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techniques for fluid problems. An algorithm for high speed flows 
was previously under development (ref. 1) when the author proposed 
low speed simulation studies in 1974. The results of the present 
investigation indicate that further research should be completed 
in order to obtain objective comparisons. The research effort of 
the ..athor was, however, redirected to enable him to participate in 
research work on "compliant walls". It is hoped that further 
extensions of the present work will be undertaken at a later date. 

THE FINITE ELEMENT METHOD 

The finite eleme it method has been used with a great deal of 
success as a simulttion technique for problems in structural 
mechanics. There are excellent books (refs. 2, 3) which detail 
the general capabilities of the method and its applications. 

However, there have been only a few known successful applications of 
this technique to fluid flow problems. Even in those cases that are 
reported in the literature no rigorous comparison has been made 
between the firite element and finite difference solutions for the 
same test problem of accuracy, efficiency, cost, or storage. It is 
felt that such a comparison should be made; hence, the test problem of 
driven cavity has been undertaken where extensive finite difference 
solutions are available (ref. 4) , 

Traditionally, the problems in fluid dynamics have been studied 
by using the "high speed variables" (primitive variables: u, v, p) 

for supersonic flows and by using the "low-speed variables" 
(streamfunction-vorticity : ¥, c) subsonic flows. In recent 

years, where mixed flows have had to be analyzed it has become 
imperative to develop numerical techniques where flow variables can 
be used at all speeds. In the literature are found the MAC, SMAC, 

ICE methods and Chorin's artificial compressibility method (refs. 5, 

6, 7, and 8) , primarily for such applications. The Fluid Mechanics 
Branch was also motivated to develop a primitive variables algorithm 
from such a consideration. 
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One of the major problems associated with incompressible flows 
is the solution of the continuity equation (hence the concept of 
stream-functions and stream tubes) . When primitive variables are 
used, the continuity equation has to be satisfied through some 
iteration procedure. A natural parameter for iteration is the 
pressure. In this report two different schemes of itei"»tion, which 
will be referred to in subsequent sections of this report as Models 
I and II, have been used. Some of the more important features of 
the present algorithm may be summarized as follows: 

• Triangular finite elements with bicubic shape functions (T-10 
element) (ref. 9), Order of accuracy 0(h**) in space. 

• Leap-frog differencing for marching to steady state (Crank- 
Nicholson scheme). Order of accuracy O(t^) in time. 

• Guyan reduction to eliminate "unnecessary” degrees of freedom. 

• A special Poisson solver for elliptic problems which conserves 
the basic conservation laws, thereby improving overall accuracy of 
the finite elements. 


Fluid Dynamic Equations 

The governing equations for two-dimensional, incompressible, 
viscous flow in a rectangular coordinate system are as follows: 
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where, u is the velocity in the x-direction, v is the velocity 
in the y-direction, and p is the pressure. 

Non-dimensionalize the flow equations using the diffusion 
scales, i.e.. 
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The non-dimensional Reynolds nmnber is defined as 
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Equations (Is, b, c) in non-dimensional form are 
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Iteration Schemes 

The system of equations (4a) to (4c) is not amenable to 
direct solution because of the nature of the incompressibility con 
dition of equation (4a) . Described below are two different models 
for the equations (4a, b, c) , 

Model I 




3u 

dT 


+ 



9u _ 9p 1 . 


32u\ 


(5a) 

(5b) 


4 


(5c) 


8v 

3i 


+ 


u 


,3v 

35 


ylZ = -lE 

^3n 3n 


Re\3C: 



where K is an artificially chosen time parameter. Equations (5a, 
b, c) describe a purely hyperbolic system when the viscous terms are 
absent and can, hence, be viewed as hyperbolic for fairly large 
Reynolds number flows (see Appendix A) . The incompressibility 
condition as modelled by Equation (5a) is a restatement, at steady 

. e . , t 

forces is equal to zero. Model I was tested on the driven cavity 
problem. It was observed that the artificial time constant "K" can 
be chosen to accelerate or decelerate the march to steady state, 

A careful choice of ”K" is essential to successfully implement 
this algorithm. At high Reynolds numbers oscillations were seen, ' 
as a direct result of ill-chosen values for "K", 
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The model equations (6) follow from the SMAC methodology (ref. 5). 

The Navier-Stokes equations are solved (6a, b) to obtain an auxiliary 
field (u*, V*); the pressure gradients appearing in equation (6a, b) 
are pseudo-pressure gradients. The pseudo-pressure is obtained by 
solving a Poisson equation (6c) . The actual velocity fields are 
obtained by applying corrections (6d, e) . The model is valid 
whenever velocity-splitting is justifiable. Assuming that this is 
so, we notice equations (6d, e) 
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u = u* - Vp 
and hence, 

7»u=V*u*-V»V*p (7a) 

By forcing 

V • u = 0 (7b) 

one obtains (6c) i.e., 

7 . Vp = V • u* (7c) 

This model was also tested on the driven cavity problem. A special 
Poisson solver, using a variational formulation (Weak-Galerkin method) 
with undetermined Lagrangian multipliers, was developed to solve the 
Poisson equation (6c) when it was found that conventional finite 
element method of satisfying the boundary conditions for this problem 
violated the basic conservation laws on the boundary. For details, 
see reference (10) . 


Finite Element Algorithm 
Triangular elements with bicubic basis functions 

o 

(J = 1, 2, . . . 10) associated with triangles are chosen in this 
work (Appendix B) . The basis functions have support only on the 
triangulation. A leap-frog differencing scheme is used on the model 
equations to yield an implicit procedure. Equation system (5a, b, c) 
after differencing can be written as follows: 
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where u*, v*, p* are the current variables, Li, L^, L 3 are the 
"residuals", the superscripts (n) and (n - 1 ) refer to previous time 
periods and commas (,) indicate partial differentials. 

It should be noted that the current variables u* , v* , p* can 

A A A 

be expressed in terms of the nodal variables u ^ , v^, p^ 

(j = 1 , 2 , . . . 10 ) as, 

u* = ? ({>.u. = {<J)}^{u} (9a) 

3 3 3 

V* = E <{».v. = {4>}'*^{v} (9b) 

j J ^ 

p* = E (J>.p. = {(|)}'^{p} (9c) 

j J 3 

The finite element algorithm is obtained by minimization of the 
residuals in equation ( 8 ) by the weights (|)j (j = 1 , 2 , , . . 10 ) 
(Galerkin -aethod) , thus yielding 30 equations in 30 unknowns 
(u, V, p) . Expansions such as 9a, b, c are used to describe the 
fields at "n" and " (n - 1 )", i.e., 

n r^,T,A.n 
P = (c|)} {p} 

(10) 

p = tp) etc. 
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After manipulations, the set of algebraic equations 


[A‘]{p> - {bi) , 

[a2]{u} - (bz) , 
[a2]{v) . {bs) . 


are obtained, where [a^] , [a^ are square matrices (10 10), 

called the element stiffness matrices, (bi), {b 2 )» (baJ "load 

vectors" (10 x l) . Details of the derivation of the stiffness 
matrix and load vector are given in Appendix C. 

Note that 


i “ ^ "" ff * "♦k,yV} 


(12a) 


i = + ■'XT * '=*k''k"> 

L J, 1 3 JJ k' 

^ S //K,X*J,X + 

n 

+ X r ( 4 . 4 . dx - dy) (12b) 

R -J, 

b'. = XT«»iPi"> - kt ffz <'fk,xV 

J JJ „ K 

+ "*k,yV} *j 

+ <4,jjV,^"> <4>i,yU."- ♦.dP. 

- I ♦j,x + ■'♦l,y"i"'^^ *j,y} 


dn 


^2 

.dn 


(12c) 




n-1. 


i 


(12d) 


= 


+ 2t ■ 2tjT <^jdy 

//<*iVi"-^> ♦jdn - ^ff^{<^k\''> <*i,x''i"'^» 

+ <Vk"* <*i,y ♦j<>“ 

- |//f*i,x''i.''‘^" ♦j.x 

- sf, hl,y''r'^ 

+ 2t //{-., p,"> .,, ydfl + 2x f <||>^P^”> 0jdx| (12e) 


It should be noted that the element stiffness ir;atrices for u, v 
have identical form, thus saving one assembly. The various integra- 
tions are done with quadrature schemes. The 16 point scheme of 
Hammer, Marlow and Stroud (11) accurate to total degree five. Is 
used for area integrals. A third order Gaussian quadrature is used 
for line integrals. To perform integrations the finite element 
triangulation is transformed to a standard triangle whose coordinates 
are (0, 0; +1, -1; +1, +1) through linear mapping T. Before any 
global assembly the stiffness matrices are transformed to the global 
values through the relation 


[Kg] = 


(13) 


Guy an Reduction 

The internal degree of freedom is often cumbersome to carry, 
incurring penalty in storage by way of increased bandwidth and 
increased number of equations. A Guyan reduction (9) is carried 
out to eliminate the internal degrees of freedom to obtain compact 
( 9 x 9 ) element stiffness matrices. 

Progrcim Documentation 

To obtain the solution of a boundary value problem, the 
following informations needs to be input to the finite element 


program: 

1. The topology of the finite element space i.e., number of 
elements « nodes, details of connectivity, flags for boundary inte- 
grals indicating the orientation of the boundary etc., 

2. The initial and boundary conditions, and 

3. Quadrature weights. 

A mesh generator was developed to create information on finite 
element topology, thereby reducing the cumbersome job of hand- 
inputting. Figure la, b shows the finite element discretisation 
obtained by using the mesh generator. Presently, the. mesh generator 
can be used only for discretisation of regular geome^i.es; however, 
extensions to irregular domains can be achieved by slight modifica- 
tions in the code. The discretisation obtained is uniform. The 
mesh generator requires a minimum number of data cards and creates 
a massive amount of data for the main program. Figure 1 was obtained 
by specifying a=l, b=l, d= 1/5, and for Figure 2 set a = 1, 
b = 1, d = 1/9 where a and b are the dim.ensions of the 
rectangular domain and d is the mesh width. 

For functional efficiency the finite element package developed 
was partitioned into two separate program segments. In Part I a new 
geometry and a new discretisation are processed and the output of 
this program consisting of quadrature values of spline function, 
matrix structure of the finite element matrices, flags for boundary 
integrals, etc. is saved. In Part II, for the given topology and 
boundary and initial conditions, finite element solutions are 
generated at each time level. All matrix operations are performed 
in Part II. The components that go into this program are 1) stiffness 
matrix assembler (to assemble element contributions into a global 
assembly) , 2) an L-U decomposition routine (where a matrix is 
decomposed into upper and lower triangular components) , 3) a matrix 
routine which converts a given matrix in local coordinate systems 
into a global matrix. 

While developing the algorithm an attempt was made to exploit 
the known spar? e structure of the finite element matrices, thereby 
enhancing savings in storage. 
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Numerical Results for the Tost Problem 


The test problem considered is a driven cavity. The upper 
wall, Figure 3, moves with uniform velocity u =* 1 for t > 0, 

Finite element solutions were obtained using Models I and II (see 
"Iteration Schemes,” pages 4 and 5). At the Reynolds numbers that 
were considered (R = 1, 10) both models behaved very well. Tables 
1 to 6 show the converged finite element fields (9x9 mesh) and 
the converged finite difference (17 x 17 mesh) solutions, for 
Re = 1.0 (Model II solutions). The finite difference solutions 
were computed using an ADI algorithm developed by Morris (ref. 12) 
Figure 4 shows the velocity profiles through the vortex center. 

Finite element approximations should show uniform convergence 
with decreasing mesh size.* To study the effect of discretization of 
solution, two finite element meshes (5 x 5, fig. 1, and 9x9, fig. 2 , 
were run. The 5x5 mesli could not capture all the details of the 
flow as was expected although the correct trend for flov; variables 
was observable. (The vorticity value at the top moving wall at 
the point midway betv/een the corners was 3.9 at Re = 1.0, rather 
low and at the point away from the singularity 9.8; again a low 
value.) The 9x9 mesh compared favorably to a 17 x 17 finite 
difference grid as can be seen from Tables 1 to 6; and was superior 
to 9 X 9 finite difference solutions. This was anticipated since 
theoretical considerations show the finite element solutions to 
have an ©(h^*) accuracy compared to the second order accurate ADI 
method. However, the finite element solutions could not reproduce 
the fourrh order accuracy. This was due to the fact that when using 
the finite element algorithm on the driven cavity, the boundary 
conditions at the singular points A and B (fig. 2) had to be modelled 
In finite difference methods the nature of differencing leaves the 
singular points alone.) With special singular function expansions, 
these singularities can be correctly modelled in finite element 
methods. The next logical step in the numerical simulation will be 


* The Tocher-10 model gives piecewise continuous first derivatives 
in the interelement boundaries and hence exemplifies a non- 
conforming element. For these elements uniform convergence as 
mesh sizes goes to zero is not always guaranteed. 
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to use the singular function expansions. Even the results that have 
been obtained through an ad-hoc modelling of the singular points 
(Tables 1-6) indicate that the finite elements can capture singulari- 
ties far better than finite difference techniques. For instance, 
comparison of the vorticities on the vertical walls (Tables 5 and 6) 
indicate the correct behavior near the singularity in the finite 
element solution that the finite difference solutions are unable to 
predict. This point is emphasized with the following logic: If 

one were to go along the moving wall towards the singularity the 
vorticity at the singularity would be +® ; along the vertical wall, 
though, the vorticity has to be negative to exhibit the driven 
nature of the flow. The corner points at the bottom should have 
zero vorticity from consideration of continuity of derivatives in 
either direction. The point on the wall closest to the singular 
point A (see fig. 3) does not therefore have the largest negative 
value of vorticity. (Since the vorticity field has to go to zero 
at the bottom point, has a negative maximum on the wall, and must 
reach +~ at the singular point, the only way this behavior can be 
exhibited by a continuous field is to have a negative maximum 
vorticity farther down). So far, to the author's knowledge no 
finite difference solution has indicated this behavior around the 
singular point. 


CONCLUSIONS 

A finite element algorithm for low speeds has been developed. 

The algorithm, applied to driven cavity problem, produced numerical 
results which were in good comparison with finite-difference results. 
Further research in this area is required in order to obtain con- 
clusive results regarding tiie efficiency of the current algorithm 
as compared to other current solution techniques such as finite 
difference procedures. 


Table 1 . Finite element solution (9x9 mesh) for velocity component 
Re. No. = 1 . 0 , converged velocity field. 
















Table 2. Finite difference solution (17 x 17 mesh) for velocity component 
(ref. 12). Converged solution. 














Table 3. Finite element solution ( 9 x 9 mesh) for velocity component 
Re. No. = 1.0, converged solution. 
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Table 4. Finite difference solution (17 x 17 mesh) for velocity component 
V (ref. 12). Re. Mo. = 1.0, converged solution. 
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Table 5. Finite element solution (9 x 9 ) of the vorticity field. Re 













Table 6. Finite difference solution (17 x 17 mesh) of the vorticity field (ref 
12); Re. No. = 1.0, converged solution. 
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Figure 2. 9x9 finite element mesh 

obtained from the mesh genera- 
tor. 
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A. PENDIJ: A 


The differential equation system modelling the incompressible 
flow is 


u, t + uu,x + uu,y ■ p,x + v(7^u) 

v, t + w,x + vv,y = p,y + vCV^y) (Al) 

p,t + p(u,x) + p(v,y) = 0 


The artificial time constant K is implicit in the time derivative 
p^. We shall take the simplified form letting the viscous terms 
be absent for analysis. Rewriting (Al) in matrix form 



(A2) 


The character of equation (A2) can be studied from the property of 
the associated matrices letting 
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it is well known that the equation system (A2) is purely hyperbolic 
if the matrix 


aA t bB , 

where a and b are arbitrary numbers has real distinct eigenvalues 
i.e., if the matrix 
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^au + bv 
0 




ap 


au + bv h 
bp 0 


(A4) 


has cistlnct eigenvalues. 

The eigenvalues of (A4) are 

Xj * au + bv 


X2 = 


au + bv ^ 1 


+ •=■ /(au + bv) 2 + 4(a^ b^) 


(A5) 


X3 = - I /(ir+ bv)2 + 4(a^ + b^)p 

For p > 0 we find that for any arbitrary a, b eigenvalues are 
real and distinct. Another consequence of equation system (A2) is 
that there need not be any boundary conditions on p; i.e., the 
pressure is continuously dependent on u and v. This is; rather 
fortuitous because the pressure equation that needs to be solved 
is a Poisson type equation which exhibits a great deal of dependeno 
on the boundary conditions. 

The presence of the artificial time multiplier k does not 
alter the conclusion made above. However the convergence rate of 
the solutions is affected drastically by the value of K. 
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APPENDIX B 


BASIS FUNCTIONS FOR T-10 ELEMENTS 



A typical finite element is shown in Figure Bl. The nodes 
1, II, III, IV have coordinates (x^, y^) ; yjT^ ^ ^^III' ^III^ 

'l ^II ^III 


and 




respectively with a 


reference x-y axis. 

The transformation. 


X = 


^l ^ 


^II ^III - 


K + 


^iii - ^ii 


y = Yt + 


^ii ^iii - 2 Yi 


(Bl) 


K + 


^III ” ^11 


If’aves the triangle with ? - n origin at node I and cr^jiiJinates 
in the transformed plane of II, III and IV as (1, -i) ; (x, +1); 
(2/3, 0) respectively. 

The basis functions (J)^, ‘J’uj be determined by 

requiring that 
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= 1 at node I 

= <|>j.j. = = 0 at II, III and IV 

(B2) 

*^1'^ = ^ 11 '^ “ '^III'^ “ '**1'’^ “ ***11'’^ “ '**111'^ ^ ^ at II and III 

♦j,,C = tj-.n = = ♦jjj.C =■ 0 at 1 

or using a total of 30 algebraic equations (The basis functions 
*^I' ^II' *^111 attached to node I such that 

= 1 = at I.) Thus convection on the line 

♦j = (1 - C) [i + 5 - - jn*] 

♦ii = <1 - 5) [? - + f\ 

♦hi = 

By rotating the cyclic order in figure B1 the "standard" triangles 
are produced with origins at II and III to obtain to The 

expression for is as follows i 

♦x = - 5) (?^ - n^) . 
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APPENDIX C 


DERIVATION OF STIFFNESS MATRIX AND LOAD VECTOR 


Elsewhere, equation (12) summarized the results of Galerkin 
minimization of the residuals (equation (8) , pages 6 and 7} . We 
develop these results in the present sections. 

Equation (8a), page 6, is reproduced below; 


Li (u*) = 1 + tu*'d^ + tv^'d - ^(D + D ) 

* X y R XX yy ' 


- 1 - Tu”d - TV% + 5(D + D ) u' 

X y R XX yy 


n-1 


+ D^p*'(2t) 


(Cl) 


Figure Cl shows a typical triangulation of a domain the boundary 

r of n is singly connected. 



Figure Cl. 

Consider a triangulation PQR of (shown hatched in above figure) . 
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The shape functions (^j have support only on this triangulation. 
We write 

u* = <fr^{u} 

V* = (C2) 

p* = ♦^{pi 

where u, v, p are the nodal variables. The Galerkin technique 
applied to equation (Cl) yields 

jTLi<J>,dSJ =0 J = 1, 2, ... 10 (C3) 

n J 

a total of 10 equations. 

The known fields at time steps (n) and (n - 1) can be expanded 
in the form similar to equation (Cl) yields 


n ,T. n „ n 

U = (() u = 

jC 

n-l ,T. n-l „ n-1 

V = ♦ ^ 


(C4) 


We also note that 


D u = D <j> u = E(}) . u . 
X x^ D/X 3 


(C5) 


Applying the integral condition equation (C3) to equation (Cl) the 
following relation is obtained; 

= // [l - T D^t + I(d^^ + Dyju"-l4,.da 

- 2TJ/D p“(^.dfi (C6) 

r\ J 


j = 1, 2, . . .10 
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I 

i 

i 


f 


The right hand side of equation system (C6) consists of known 
functions and hence gives rise to the load vector. The left hand 
side of (C6) is an algebraic expression in the nodal variables u. 

The equation system (C6) can be summarized as 

EAijUi = bj , j = 1, 2 ... 10 (C7) 

where 



= ff + T <Vk“^ 

n 

-5 (Is *(ly *i,y)*j] 


(C8) 


Due to the nature of the finite element that we chose, (piecewise 
continuous first derivatives across interelement boundaries) we 
apply Green's theorem to modify the last two terms in right hand 
side of equation (C8) , which is stated below. 


yy*VxFdfl = y Ends (C9) 

n r 


Thus, equation (C8) modifies to 


A. . 

X 1 




The weak form of Galerkin minimization procedure enables us to 

order the element connectivity such that the boundary integrals 

cancel out for all interior nodes. Only for boundary nodes need the 

boundary integrations be carried out. Similar relation as (CIO) can 

be found for b . . 
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